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Abstract 

We study a class of swarming problems wherein particles evolve dynamically via 
pairwise interaction potentials and a velocity selection mechanism. We find that 
the swarming system undergoes various changes of state as a function of the self- 
propulsion and interaction potential parameters. In this paper, we utilize a pro- 
cedure which, in a definitive way, connects a class of individual-based models to 
their continuum formulations and determine criteria for the validity of the latter. 
H-stability of the interaction potential plays a fundamental role in determining both 
the validity of the continuum approximation and the nature of the aggregation state 
transitions. We perform a linear stability analysis of the continuum model and com- 
pare the results to the simulations of the individual-based one. 
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1 Introduction 



The collective behaviors of aggregating organisms are of interest in various 
fields, including biology, engineering, mathematics, and physics [1,2,3,4,5,6,7,8,9] 
There are primarily two classes of pertinent models: individual-based and con- 
tinuum ones. In the first case, one considers a collection of individual enti- 
ties, so that the system is defined on the "microscopic" scale. Such models are 
particularly useful for the study and algorithmic design of small-size aggregates 
such as artificial swarms of autonomous vehicles. Larger discrete systems are 
adaptive to statistical analysis [1,10,11,12,13,14,15,16,17,18,19,20,21]. Contin- 
uum models typically describe swarms through a density function p (r) and a 
velocity vector field u{f). These obey appropriate non-linear and often non- 
local equations. One may presume these equations are derived from, or at least 
connected back to, the original microscopic system. Continuum models are 
useful for theoretical analysis of swarming systems [2,3,14,15,22,23,24,25,26]. 
Although both individual-based and continuum models are applicable, the 
connection between the two has, in fact, not been particularly well estab- 
lished. A primary purpose of this paper is to better unify the two approaches, 
for a particular class of models, following the classical statistical studies of 
fluids [27]. We investigate the validity of the continuum model by a detailed 
comparison with the associated individual based one. In particular, for certain 
interaction forms, the two descriptions yield the same morphological patterns. 
Furthermore, and perhaps more importantly, we are able to explain why the 
continuum model fails qualitatively for the other cases, where discrepancies 
exist. 

In Rcf. [28], a criterion from classical statistical mechanics known as H-stahility 
was applied to individual based swarming models. A system of interact- 
ing particles is said to be H-stable if the potential energy per particle is 
bounded below by a constant which is independent of the number of par- 
ticles present [29]. H-stability is a necessary and sufficient condition for the 
existence of thermodynamics. Indeed, a system without this stability will, in 
the thermodynamic limit, collapse onto itself; such systems are called catas- 
trophic. In Ref. [28], numerical simulations strongly suggest that a specific 
non-Hamiltonian swarming system exhibits the same stability trends observed 
in classical Hamiltonian many-body systems. In this paper, we show that H- 
stability also plays an important role in determining the correct passage to 
the continuum limit. 

This paper is organized as follows. In Sec. 2, the individual-based model is pre- 
sented. We study various aggregation states and transitions between them via 
numerical simulations. In Sec. 3, a continuum model is derived. In Sec. 4, we 
quantitatively compare steady states of both continuum and discrete models. 
We show that, while the proposed continuum model works well in the catas- 
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Fig. 1. Left: The swarming pattern of a single mill. Right: The swarming pattern of 
two interlocking mills. 

trophic regime, discrepancies arise for large H-stable systems. In Sec. 5, the 
stability of the homogeneous solution of the continuum model is studied and 
compared to the numerical results of the individual-based model. In Sec. 6, we 
discuss the choices between soft-core and hard-core interaction potentials. 



2 The individual-based model 

Background 

Common swarming patterns have been observed and reported in various species 
in nature. One example is a coherent flock formation involving a polarized 
group moving in the same direction. Another example is a single rotating 
mill pattern, with a rather stationary center of mass, as in Fig. 1 (left). The 
rotating-mill pattern is frequently observed in both two and three dimensions 
among many species and across different sizes [5,7,30]. Various individual- 
based models have been able to reproduce these patterns within certain pa- 
rameter ranges [13,15,16,17,18,21]. An unusual pattern of overlapping double 
mills is also reported in Ref. [15], similar to the simulation shown in Fig. 1 
(right). The double-mill phenomenon is observed in the early stages of aggre- 
gation of Myxococcus xanthus, a single-cell bacteria driven by self-propelling 
motors [31]. We show that this configuration can be obtained using the same 
swarming mechanism that produces the single-mill pattern but exists in a dif- 
ferent parameter regime. The rarity of the double-mill state is discussed in 
Sec. 6. 
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Equations of motion 



The swarming model we present in this paper is described by the following 
equations of motion 



mi-^ = avi- P\vifvi-VUi, (2) 



where rrii, Xi and Vi are, respectively, the mass, position, and velocity of particle 
i. The terms avi and — /5 Iwil^Wj define the mechanism of self-acceleration and 
deceleration which give the particles a tendency to approach an equilibrium 
speed Veq — \joi/ (5. This Rayleigh-type dissipation was originally proposed in 
Ref. [32] and is often used in the literature as a velocity-selecting mechanism 
[2,10,13,17,18,21,33]. The potential Ui describes the interaction of particle i 
with the other particles. One common choice is the following [15,18,23,28] 



Ui = U {xi)=Y.V (If, - = ( -Cae-^^ +C',e-^^ J . (3) 

Eq. (3) assumes that only pairwise interactions are significant and ignores N- 
body interactions with > 3. The pairwise interaction consists of an attrac- 
tion and a repulsion with Ca, C,. specifying their respective strengths and 
Ir their effective interaction length scales. Similar behaviors are also observed 
with other functional forms of interaction potential that are characteristically 
similar to Eq. (3). Note that to simplify the analysis, our model is determin- 
istic. Stochastic forces appear in many other models [1,10,13,17,18,21]. In our 
simulations, we observe that noise affects the swarming patterns only beyond 
certain thresholds, and thus its consequences are not investigated. 

We can non-dimensionalize the equations of motion by substituting t' — 
(rni/£a^(3^ t, x\ = f and thus, = {ia^jm^ Vi into Eqs. (1) - (3) 



f = «'^-|^r^-;^V,^', (5) 

U/^J2(-^~^''~'''^ +Ce-^Y (6) 
j^i \ J 

where a' = (a/34') /m^, mi = m^V (fCJa). C = C./a, and £ = 4/4; 
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Fig. 2. The H-stability diagram of the interaction potential in Eq. (3) [28]. The 
shaded region is the so-cahed biologically relevant region where the interaction con- 
sists of a long-range attraction and a short-range repulsion. 



hence, the model is essentially a 4-parameter one. In Ref. [28] the effects of 
varying C and ^, which affect H-stability, are explored. In this paper we inves- 
tigate the role of a', the relative strength of the self-driving force with respect 
to the interaction. The parameter m/ affects the time scale of the particle 
interaction and is fixed during our investigations. Note that the dimensional 
parameter a only appears in the dimensionless parameter a', which allows us 
to vary a to change a' without affecting the other three independent parame- 
ters, provided that /3, ^a, and rrii are fixed during the process. To preserve the 
original meaning of the model parameters, our results are presented in the di- 
mensional form by using Eqs. (1) - (3). Only the biologically relevant cases that 
consist of a long-range attraction and a short-range repulsion are studied. In 
other words, we confine our analysis to the parameter space where C > 1 and 
^ < 1, which is shown by the shaded region of the H-stability phase diagram 
in Fig. 2. The extremely collapsing cases reported in [28], such as the ring for- 
mations and the clump formations illustrated in other regions, do not change 
morphology with respect to a' . 
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Swarming states 

We use the fourth order Runge-Kutta and the four step Adam-Bashforth meth- 
ods for the numerical simulation of Eqs. (1) - (3) [34] . We impose free boundary 
conditions to the model and initiate the simulation with random distributions 
of particle position and velocity. Figure 1 shows two typical patterns akin to 
those observed in various natural swarms. On the left panel is the single-mill 
state, where every particle travels at the same speed Veq around an empty 
core at the center of the swarm. On the right panel is the double-mill state, in 
which particles travel in both clockwise and counterclockwise directions, also 
at a uniform speed fcq. In this second example, when viewed as two superim- 
posed mills, the cores of each mill do not exactly coincide but rather fluctuate 
near each other. Another two states are shown in Fig. 3. On the left panel is 
the coherent flock state. All particles travel at a unified velocity (i.e., with the 
same speed and direction) while self-organizing into a stable formation. On the 
right panel is the rigid-body rotation state. The flock formation closely resem- 
bles that of the coherent flock, but instead of traveling at the same velocity, 
the particles circulate around the swarm center defining a constant angular 
velocity co. Unhke the single and double-mill state, where particles swim freely 
within the swarm, both the coherent flock and the rigid-body rotation states 
bind particles at flxed relative positions, exhibiting a lattice-type formation. 
Hence, we also use the term lattice states to refer to both the coherent flock 
and the rigid-body rotation states. Note that the coherent flock is a traveling 
wave solution of the model, and thus a solution of the following Euler-Lagrange 
equation 

W, = V^V(|f,-f,|) = 0. 

It is interesting to note that this equation arises in the context of a gradient 
flow algorithm for autonomous vehicle control [8,35,36,37]. Thus the flock 
formations have the shape and structure as equilibria of the gradient flow 
problem with the same potential. 

The coherent flock and the single-mill states are among the most common 
patterns observed in biological swarms [5,7,30]. The double- mill pattern is 
also occasionally seen; an example is the M. xanthus bacteria at the onset of 
fruiting body formation [31]. On the other hand, natural occurrences of rigid- 
body rotation, to the best of our knowledge, have not been reported in the 
literature. Indeed, the rigid-body rotation, where every particle travels at a 
constant angular velocity uo, does not deflne a rotationally symmetric solution 
for Eqs. (1)- (3) and the swarm is observed to drift randomly due to the un- 
balanced self-driving mechanism. The random drift may eventually break the 
rotational symmetry and turn the swarm into a coherent flock after a tran- 
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Fig. 3. Left: The coherent flock state. Right: The rigid-body rotation state. 

sient period. Thus, we speculate that this pattern may only be a meta-stable 
or a transient state. In addition to the above aggregation states, the particles 
may simply escape from the collective potential field, and no aggregation is 
observed. We name it the dispersed state. 

Using numerical simulations, we find that the H-stable swarms undergo a 
different state transition process from that of the non-H-stable swarms. For 
both H-stable and catastrophic interactions, the lattice states, as shown in 
Fig. 3, emerge for low values of a, and thus, of low t>eq. In this case, the 
confining interaction potential is stronger than the kinetic energy of individual 
particles and tends to bind the particles at specific "crystal" lattice sites. Most 
initial conditions lead to the coherent fiock state while some occasions result 
in the rigid-body rotation state. The state transition of H-stable swarms is 
simpler. As a increases, the particles eventually gain enough kinetic energy to 
dissolve the aggregation. 

The state transition of catastrophic swarms is characterized by more behav- 
ioral stages. Starting from the lattice states and upon increasing a, the parti- 
cles gain more kinetic energy from the environment to reach f eq and are able 
to break away from the crystal lattice sites. However, unlike H-stable swarms, 
the interaction potential in the catastrophic regime is still strong enough to 
aggregate medium-speed particles within a swarm. In this regime, core-free 
mill states emerge, as shown in Fig. 1. Since all particles travel at a non-zero 
uniform speed, the centripetal force provided by the collective interaction po- 
tential is not strong enough to sustain such particles too close to the rotational 
center. As a result, the mill core is a particle-free region. At moderate a, a 
single mill state emerges. At slightly higher a, we observe both single mills 
and double mills as possible states. In the latter case, the interaction potential 
gradually loses its effectiveness to unify the clockwise (CW) and counterclock- 
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wise (CCW) rotational directions; particles traveling in the opposite direction 
with respect to the majority tend to not change their direction of motion, 
and double mills can emerge. The transition from single to double mill is a 
gradual process. Figure 4 shows the number of particles in each rotational di- 
rection for various values of a. In the single-mill regime, particles traveling 
at one direction are quickly assimilated into the other (Fig. 4, top). Upon in- 
creasing a, the particles no longer settle into a unified rotational direction 
(Fig. 4, middle), and for large enough a, approximately the same number of 
particles travel in each of the CW and CCW directions (Fig. 4, bottom). The 
presence of either a velocity alignment rule or a hard-core repulsive interac- 
tion will destroy this double-mill state. The latter case is because hard-cores 
always provide a system with H-stability. Thus, it is clear that for sufficiently 
many particles, the double mills will ultimately break apart. Notwithstanding, 
it appears that the double mills are especially sensitive to hard-cores and, even 
for small cores and moderate N , we have not observed these structures. As for 
the coherent flock state, it still remains a possibility in this regime where the 
mill states occur. However, the basin of attraction is greatly reduced, and only 
very polarized initial conditions can lead to the coherent flock formation. As 
a increases beyond the double-mill regime, particle kinetic energy eventually 
becomes high enough to break up the swarm. This is the dispersed state, and 
no aggregation can be found. 

Upon fixing the other parameters, the threshold between the aggregation and 
the dispersed states is described by a critical escape value of a, denoted by ctesc- 
Figure 5 shows a^sc of an H-stable swarm versus a catastrophic one, in which 
single-mill states are generated and a is increased until the dispersed regime is 
attained. For the H-stable case, Ocsc does not vary significantly with respect to 
the total particle number of the swarm, denoted by iV, due to the fact that the 
nearest neighbor distance (5nnd) does not change as N increases. As a result, 
the binding potential energy of the interaction force acting over each particle 
is independent of N . On the other hand, ctesc of the catastrophic swarm varies 
linearly with respect to N . From our numerical simulations, we observe that 
the outer and the inner radii of the catastrophic swarm remain approximately 
fixed with respect to N while a < a^sc- Based on this observation, we can 
derive a semi-empirical formula to estimate the value of a^sc by assuming that 
particles are uniformly distributed in a doughnut shape domain. By balancing 
the centripetal and the interaction forces, we obtain 



A T Rout 
Win 

where i?in and -Rout denote the inner and the outer radii of the single mill, 
respectively, and x is an arbitrary unit vector. This estimate predicts that ctesc 
should scale linearly with N, which is clearly illustrated in Fig. 5, where we 
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Fig. 4. Time variation of the numbers of particles rotating in different directions: 
Tlie triangles represent the number of CCW particles while the circles are of CW 
ones. (Top) a = 1.5. (Middle) a = 4.0. (Lower) a = 6.0. The fixed parameters are 
P = 0.5, Ca = 0.5, Cr = 1.0, ia = 2.0, ir = 0.5, and N = 500. All parameters here 
and throughout the paper are in arbitrary units. 

use the numerically simulated Rout = 5.2 and Ria = 1.2 for a quantitative 
comparison. 



State transitions of H-stable and catastrophic swarms 

In order to quantitatively determine whether the swarm is in a coherent flock 
state or a single-mill state, Couzin et al. have proposed two measures [16]: the 
polarity, P, and the (normalized) angular momentum, M, defined as follows 
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Fig. 5. Oesc versus the total number of particles in an H-stable swarm (/3 = 0.5, 
Ca = 0.5, Cr = 1.0, £a = 2.0, ir = 1-5, dashed line) compared to that of a catas- 
trophic swarm (/? = 0.5, Ca = 0.5, Cr = 1.0, 4 = 2.0, 4 = 0.5, dotted line). The 
solid line is the curve estimated by Eq. (7). 
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(9) 



where rj = Xj — xcm, and xcm is the position of the center of mass. A perfect 
coherent flock results in P = 1 and M = while a perfect single-mill pattern 
results in M = 1 and P = 0. In order to distinguish the double-mill pat- 
tern, we propose an additional measure by modifying the normalized angular 
momentum 
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(10) 



If a double-mill pattern has perfectly equal numbers of particles going at each 
direction with the centers of mass of both directions exactly overlap, Mabs = 1 
and M = 0; both M and Mabs equal to 1 for a single mill. 
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Fig. 6. Emergence of a rotating single-mill structure from a rigid-body rotation 
in the catastrophic regime. The left panel shows the ensemble averaged tangen- 
tial velocity, {v (r))tang, of particles at a distance r from the center of mass. Each 
('^))tang figure corresponds to the swarm structure of different values of a on the 
right panel: from top to bottom are a = 0.003, a = 0.03, a = 0.1, and a = 0.5. The 
other parameters are /? = 0.5, Ca = 0.5, C,. = 1.0, 4 = 2.0, 4 = 0.5, and N = 500. 



11 



Although the presence of the coherent flock that yields P ~ 1 allows us to use 
P to quantify the transition from lattice to single-mill state, the co-existing 
rigid-body rotation state, for which P ~ 0, introduces spurious events. Since 
the rigid body state has a much smaller basin of attraction than the coherent 
flock, one choice is discarding all rigid-body rotation events and selecting only 
the coherent flock ones. However, the boundary between a rigid body rotation 
and a single mill is ambiguous, as shown in Fig. 6, where a rigid-body rotation 
transforms to a single mill by increasing a. Since a constant tangential speed 
indicates a milling formation, and a constant angular velocity (i.e., a linear 
tangential speed against r) characterizes a rigid-body rotation, we can see 
from the figure that two states are mixed during the transition: the outer 
part of the swarm begins to exhibit the milling phenomena while the inner 
part still remains a rigid body. Indeed, the collective interaction potential is 
stronger in the inner part of the swarm, and particles need a higher kinetic 
energy injection from the self-driving terms to escape the binding potential. 
Since the lattice formation of the rigid-body rotation has an ordered particle 
distribution, and the milling formation exhibits more disordered distribution, 
we propose an ordering factor of period Q to quantitatively distinguish these 
two states 



N M 



=1 j 



(11) 



where is the angle between Xij and iCjj+i with defined as Xj — Xi. 

The summation index j here represents the j-th nearest neighbor of particle 
i, and n denotes the number of neighbors that are taken into consideration for 
each particle. We also define Xi^^+i = Xi^i to simplify the formula. If all 4>fj+i 
are distributed at 2nk/Q where < Q is a positive integer, 0(q) = 1, and 
the particles are distributed on a lattice of period Q. On the other hand, if 
the distribution is completely random, cancellation occurs in the summation 
of cosines, and ~ for all Q. The number of nearest neighbors of each 
particle i can be arbitrarily chosen for // > 2. However, note that cannot 
be too large; otherwise, second layer neighbors may be counted, which results 
in an incorrect Q. For the sake of definiteness, we choose /i = 3. In order to 
avoid incorrect estimations due to the dispersed state, we also impose that 
a particle pair must be separated by a distance no larger than 2£a for the 
particles to qualify as neighbors. Figure 7 (b) shows the distribution of 0^ j+i 
collected for all i and j on a rigid-body formation. Peaks are observed at kn/S 
{1 < k < 5), indicating that the formation is a hexagonal lattice. Figure 7(c) 
shows 0(Q) versus Q for the same rigid-body formation. As expected for a 
hexagonal lattice, the curve peaks at Q = 6. Therefore, 0(6) can be used to 
explore the transition from a hexagonal lattice to a non-lattice mill state. 
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Fig. 7. (a) The ordering factor of period 6 versus a and an illustration showing the 
definition of (/'j j^^. The triangles are data points of the catastrophic case while the 
squares represent the H-stable case. The parameters other than q for both cases 
are the same as those in Fig. 5 with = 200. (b) The distribution of for all 

i and j. (c) Comparison of the ordering factors of different periods Q. 



Using the quantities defined in Eqs. (8) - (11), different swarming states can be 
classified. Dramatic changes in P, M, Mabs, and 0(q) are observed upon mod- 
ifying specific parameters in the model and indicate a change in the swarming 
state. Figure 7 (a) shows the transition from lattice to single-mill states for 
catastrophic swarms as 0(6) gradually decreases with respect to increasing a. 
Also shown in the figure are the same quantities for an H-stable swarm; note 
that as a increases, 0(6) suddenly drops to zero, corresponding to the sudden 
dissolution of the hexagonal lattice structure into a dispersed state. The larger 
value of 0(6) in the H-stable swarm indicates a more regular hexagonal lattice 
formation. 



For higher values of a, we further consider P, M, and M^bs to differentiate the 
coherent state and the two mill states. Additionally, in order to distinguish 
the dispersed state from the rest, we calculate the aggregation fraction, /agg, 
defined as the fraction of the initial particles that aggregate as a swarm. 
In Fig. 8, we show how H-stable swarms differ from catastrophic ones during 
the transition between states. Figure 8 (a) shows that the H-stable swarm is a 
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Fig. 8. The state transition diagram of (a) an H-stable swarm and (b) a catastrophic 
swarm. The fixed parameters are the same as those in Fig. 5 with = 200. 

coherent flock for small a, indicated by P ~ 1. For increasing a, the swarm 
disperses and /agg = 0. Note that P remains close to one when /agg ^ 0, 
indicating that the aggregate goes from the coherent lattice state directly to 
the dispersed one. Figure 8 (b) shows the transition of a catastrophic swarm, 
which displays a full four-stage transition: in the small a regime, particles 
arrange as a coherent lattice with P ~ 1; as a keeps increasing, the single-mill 
state appears (P ~ and M ^ 1), followed by the double- mill state (Mabs — 1 
and M ~ 0) until the dispersed state (/agg = 0) is reached. 

Drawing an analogy from the state transition of swarming patterns to the 
phase transition of materials, the lattice states can be regarded as "solid" since 
interparticle distances are kept constant. The milling state allows particles to 
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"swim" within a finite volume without being bound to a fixed lattice site; thus, 
it can be regarded as "liquid" . Finally, in the dispersed state, particles escape, 
similarly to a "gas". Upon increasing a, a catastrophic swarm undergoes the 
solid-liquid and liquid-gas transitions, which resemble the processes of melting 
and vaporization. On the other hand, an H-stable swarm goes from the solid 
state directly to the gas one, which is more similar to sublimation. 

Consistently with granular media models [38], we may define a "temperature" 
analog using the variation of the individual particle velocity among the flock: 
Tg = (^{v — {v))'^Y Note that (v) is the velocity of the center of mass, and 
thus, Ts ~ for the coherent flock pattern, while Tg ^ for the steady mill 
states. The swarming patterns change from one state to another by varying 
T 

While different aggregation morphologies can be studied using the individual- 
based model, the large number of degrees of freedom involved pose a difficulty 
for analyzing the dynamics of large systems. In the following section, we 
therefore develop and investigate a continuum model consistent with the mi- 
croscopic description of Eqs. (1) - (3). 



3 Continuum model 



The continuum approach is widely adopted for modeling swarming systems, 
especially on the ecological scale where massive movements of populations 
are considered [2,3,22,23,25,26]. It is also more suitable for theoretical anal- 
ysis especially in large N limit. Due to the lack of connections between in- 
dividual rules and continuum "fluxes" , most continuum models in the litera- 
ture are constructed on the basis of heuristic arguments. Attempts have been 
made to bridge the gap. For example in Ref . [39] , a continuum kinematic one- 
dimensional advection-diffusion model is derived based on a biased random 
walk process of a set of Poisson-distributed particles. This work is extended to 
higher dimensions and to more general kinematic rules in Ref. [40,41]. However, 
our model is based on dynamic rules and the corresponding continuum limit 
is much more difficult to rigorously justify. In Ref. [14] , a continuum model 
is derived from a class of dynamic individual-based descriptions by using a 
Fokker-Planck approach. In order for the fiux term to be amenable to analytic 
investigations, the Fokker-Planck equations have to be closed under several 
assumptions, but these assumptions, such as that the preferred velocity which 
particles tend to reach is small with respect to noise terms, are not applicable 
to our model. 

In this paper, we derive a continuum model by explicitly calculating the en- 
semble average the model of Eqs. (1) - (2) using a probability distribution func- 
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tion. This classical procedure is described in Ref. [27] where continuum hy- 
drodynamics equations are derived starting from a microscopic collection of 
particles. Let 

/ ^ f(xi,X2,...,XN;Pi,P2,-,PN;t) (12) 



be the probability distribution function on the phase space, defined by position 
and momentum [xi^pi), 1 < i < N, at time t. The mass density p{x,t), the 
ensemble velocity field u{x,t), and the continuum interaction force Fy {x,t) 
can be defined as 



N 



p{x,t)^m^{5{xi-x)]f), (13) 



i=l 



^ ' Pix,t) pix,t) ^ ' 

N 

Fv (f , i) = E (-Vx- C/ {xi) 5 {xi - x) ■ /). (15) 

1=1 

We consider the case of identical masses, mi = m. The function 6 (x) is the 
Dirac delta function, and U (x,) the collective interaction potential acting 
on particle i. Using the generalized Liouville theorem that incorporates the 
deformation of phase space due to the non-Hamiltonian nature of the system 
at hand [42], we obtain the continuum equations of motion 



| + V.(p.l^O, 



(16) 



d_ 
di 



(pu) + V • (puu) + V -aK^apu- 2(3Eku - 2/3gK + 2/3w • gk + Fy. (17) 



The first is the equation of continuity, and the second is the momentum trans- 
port equation. Here, E-^ = p\uf /2 is the kinetic energy. The terms qK{x,t) 
and (3"k {x, t) are mathematically defined as 



N 



qk {x, t)^J2 



m 



1=1 

N 



Pi 



m 



aK{x,t)^^mi [ — -u \ [ — - u\5{xi- x)]f), 



u 



Pi_ 

m 

'Pi 



u\5{xi-x)-j), 



i=l 



m 



m 



and represent the energy flux and the stress tensor due to local fluctuations in 
particle velocities with respect to u{x, t). The derivation of the term V-(5"k can 
be found in Ref. [27] while the other terms related to and ax are derived in 
Appendix A. By simulating the discrete model, we estimate the magnitude of 
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gk and o"k and find that both fluctuation terms become negligible with respect 
to the the other terms on the RHS of Eq. (17) in the lattice, single-mill, and 
the dispersed states. Thus, neglecting the fluctuation terms, we obtain 



|^ + V-W=0, (18) 

— (pu) + V • (puu) = apu - 2/3Eku + Fy. (19) 
at 



Continuum interaction force 



In Eq. (19), the continuum interaction force can be obtained by substituting 
the exphcit form of the interaction potential Eq. (3) into Eq. (15) 

N N 

Fv (f , t) = E E {-^^y - ^j) ^ - ^) ; /)■ (20) 

1=1 j=l 



Using the fact that an arbitrary function F (xj) Vx, e R'' can be written as 
F{xj)^ J F{y)S{xj-y)dy, 

we can rewrite Eq. (20) as 

N N „ 

Fy(f,t)=EE / dy{-Vg^V{x,-y)6{x,-x)6{x,-y)-f) 

„ N N 

= / -V,V {x-y)Y.Y.{^ - V) S {x, - x) ; f)dy 

= J -VgV {x - y) (a;, y, t) dy, (21) 

where the p^"^^ is the pair density 



N N 

(f , ^, t) ^ ^ ^ (5 {x, - y) 6 {x, -x)-f). 

i=i j=i 



Note that we should take the ensemble average on a scale considerably larger 
than the spacing between particles. If the particles are quite dispersed, the 
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suitable scale may be much larger than the characteristic lengths of the in- 
teraction force {—'VV in Eq. (21)), rendering it localized. In this case, the 
continuum approach cannot capture the swarming characteristics occurring 
on the the interaction scale and fails to describe the individual-based model 
on such a scale. This is what occurs in the H-stable regime, which we further 
discuss in Sec. 4. 

For identical particles, the pair density can be written as 



m 



p(^) (f , y, t) = -^p (f , t) p {y, t) ^(^^ {x, y) , 



where the correlation function g^"^^ {x, y) — 1 when the particles have no in- 
trinsic correlation. Using this assumption, 

p(2)(f,y,i) = ^p(f,i)p(y,i), (22) 



and 



Fv (f , t)= j -W^V {x - y) \p (f , t) p {y, t)dy. (23) 



If we further substitute the interaction potential specified in Eq. (3) into the 
above equation, we get 



Fy {x, t) = -p {x, t)V j (-^ e^^ +^ e-^) p {y, t)dy. (24) 



Since we assume that all particles have an identical mass, we may choose m = 1 
without loss of generality. In this case, Eq. (24) becomes the one proposed in 
Ref. [15], on heuristic grounds. Using Eq. (18), we may modify Eq. (19) and 
divide by p on both sides to obtain a more conventional expression 



^ + V-{pu)^0, (25) 
+ u ■ Vu — au — (3 \uf' u ( V {x — y) p {y, t)dy. (26) 

C/b TJX J 
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4 Compcirison to the individual-based model 



The time-dependent variations of the density p (x, t) and of the momentum 
p (5, t) = p (x, t) u (x, t) can be obtained through numerical simulations of 
Eqs. (18) -(19). Here, we use the Lax-Friedrichs method [43] to integrate the 
partial differential equations. We compare the results of the continuum model 
to those of the individual-based model of Eqs. (l)-(2). Figure 9 shows the 
frequently observed single-mill steady state solutions of both models in the 
catastrophic regime. Both simulations use identical parameter values and the 
same total mass 



which are listed in the figure caption. In the individual-based model, particles 
are initially distributed with random velocities and at random positions in a 
2£a X 2£a box. The initial condition of the continuum model is a homogeneous 
density in a box of the same size with randomized momentum field. While 
the individual-based model adopts a free boundary condition that allows the 
particles to move around over the entire space, the continuum model uses 
an equivalent boundary condition of out-going waves but on a fixed compu- 
tational domain of 5£a x 5£a size which is divided into 256 x 256 grid cells. 
The individual-based simulation uses an adaptive time step size that keeps 
the increment in position of each step under £a/10 and the increment in ve- 
locity under Ca/5m. The time step size of the continuum simulation is chosen 
so that the CFL number is 0.98. Figure 9 (a) illustrates the averaged den- 
sity (p) as a function of the radial distance from the center of mass. These 
two profiles are in good agreement despite the density oscillation shown in 
the individual-based model, reflecting a multiple-ring ordering of the particle 
distribution. Figures 9 (b) and (c) match the averaged radial and tangential 
momenta (denoted by (p)rad and (p)tang respectively) from the simulations of 
both models. The negligible radial momenta in Fig. 9 (b) indicate that there is 
no net inward or outward mass movement, and thus, the density profile along 
the radial direction is steady. We can divide the momentum by the density 
to obtain the velocity field. In Fig. 9 (d), we show the averaged tangential ve- 
locities, (w)tang = (p)tang/(p), from the simulations of both models; it shows 
that both the individual-based and the continuum swarms are rotating at the 
same constant speed, which equals to v^. 
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Fig. 9. Comparison of the numerical simulations of the individual-based model and 
the continuum model: The parameters used in both simulations are Ca = 0.5, 
Cr = 1.0, ia = 2.0, £r = 0.5, a = 1.2, P = 0.5, and the total mass mtot = 88. (a) 
The averaged density profiles along the radial distance from the center of mass, (b) 
The averaged radial momentum profiles, (c) The averaged tangential momentum 
profiles, (d) The averaged tangential velocity profiles. 



Validity of the continuum model 



The ensemble average implicit in the continuum approach does not allow for 
double-milling in the continuum limit because the velocity inside a mesh cell 
is averaged and unified. Local velocity variations, which contribute to gx {x, t) 
and aK{x,t) in Eq. (17), are neglected. We calculate the ratio of the speed 
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Fig. 10. Relative speed fluctuations Ak while forming a single mill from random 
initial conditions. The dashed curve illustrates the normalized angular momentum 
M in Eq. (9) while the solid curve represents Ak- The parameters of this simulation 
are a = 1.0, /3 = 0.5, Ca = 0.5, Cr = 1.0, £a = 2.0, = 0.5, and N = 500. 



variation to the equilibrium speed, Ak = y {{v — Veq) )/veq^, in order to ef- 
ficiently estimate the contribution of these local velocity fluctuation terms. 
Figure 10 shows that Ak becomes negligible after the swarm has reached the 
single-mill configuration and M ^ 1. However, during the transient time, Ak 
is significantly larger, which implies that {x, t) and (Jk {x, t) cannot be ne- 
glected during this period. Hence, the continuum model of Eqs. (18) - (19) can 
be useful in analyzing the stability of the steady state solution but does not 
capture the dynamics of the swarm settling into this steady state. 

While Figure 9 shows good agreement between the steady state solutions of 
the continuum and the individual-based models in the catastrophic regime, 
inconsistencies arise as the parameters shift into the H-stable regime. Here, 
at low particle speeds, the individual-based model results in compactly sup- 
ported solutions similar to those shown in Fig. 3. Conversely, the corresponding 
continuum model yields a uniform density distribution spreading over the en- 
tire computational domain regardless of domain size. Since a valid continuum 
model should reflect the large number limit of the individual-based model, we 
can investigate how the steady state solutions of the individual-based model 
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(a) 





Fig. 11. (a) Flock radius versus number of particles with the total mass fixed at 
rritot = = 500. Solid circles represent the H-stable flock {£r = 1-5), fitted by 
R oc A^*^'^^, while the solid diamonds represent the catastrophic flock {ir = 1-3), 
fitted by i? oc A^*^'^^. The other parameters are a = 0, f3 = 0.5, Ca = 0.5, Cr = 1.0, 
and ia = 2.0. (b) Exponents Z of the power law fitting R oc for a range of C 
and i. The dimensionless parameters C and £ are changed by varying Cr and £r 
while the other parameters remain the same as above. The color map indicates the 
power Z, and darker colors represent higher exponents. The solid curve marks the 
boundary between the H-stable and the catastrophic regions. 
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evolve by increasing A'^ while keeping all continuum variables and parameters 
fixed. In particular, we increase while keeping the macroscopic parameter 
^tot = fixed. If the individual-based solution has converged to the contin- 
uum limit, the solutions should be independent of any microscopic parameter, 
such as N. In Fig. 11 (a), we show the radius R of steady swarms versus N 
under fixed mtot for an H-stable case and for a catastrophic one. The flock 
size is indeed independent of for catastrophic swarms, and the two models 
yield consistent solutions, as shown in Fig. 9. However, in the H-stable regime, 
the swarm size increases with y/N in spite of a fixed mtot- This suggests that 
a compactly supported solution does not exist in the large number limit of 
an H-stable swarm. Figure 11(b) further illustrates this point by expanding 
the investigation to a broader parameter space. The H-stability threshold is 
the solid curve, parting the C-i phase space in Fig. 11 (b). The flock radius 
R is approximately independent of N for catastrophic swarms, but when the 
parameters C and i cross over to the H-stable regime, R scales as with 
Z~l/2. AsA/"— >oo, H-stable swarms tend to occupy the entire space. 

The cue to the inconsistency between the solutions of the two models in the 
H-stable regime lies in the derivation of the continuum model. As previously 
mentioned, the macroscopic variables are obtained as ensemble averages over 
a large number of microscopic ones. In the catastrophic regime, 5nnd ^ ^a, C-r 
in large A^ limit, as shown in Fig. 11. Hence, as A^ ^ oo, the particle distribu- 
tion converges to a continuum density on a scale comparable to the interaction 
length. On the other hand, for an H-stable swarm, 5nnd stays non-negligible 
with respect to the characteristic length of the interaction. Hence, Eq. (21) 
does not hold on a scale comparable to the interaction length, and 
suit, Eq. (23) is not a valid description of the continuum force on such a scale 
in the H-stable regime. This is also verifled in Fig. 12. In Fig. 12 (a), we de- 
fine significant neighbors of a particle as neighbors that exhibit a "significant 
interaction". The quantitative definition is illustrated by the graph on the 
upper-right corner, in which the pairwise interaction potential V (r) is plotted 
versus the inter-particle distance r, and the potential well depth is denoted by 
Knin- We define a distance so that V {x) > sV^in if a; > Xg, where < s < 1 
is a ratio. Then, the number of significant neighbors of each particle is the 
number of neighbors at a distance x for which x < Xg- In Fig. 12 (a), we count 
the averaged number of significant neighbors of a particle, denoted by n^, for 
s = 0.5 and s = 0.1. In the H-stable regime, Hs is very low and remains steady; 
it rises rapidly when the parameter crosses over into the catastrophic regime. 
The results suggest that the H-stable swarms are locally too sparse for Eq. (21) 
(and thus, Eq. (23)) to remain valid. Furthermore, we can use ensemble aver- 
ages to approximate the collective interaction potentials in the two models. If 
the continuum limit properly describe the individual-based description, these 
two potential energies should converge as A^ increases. Let us define Ueu as 
the continuum ensemble average interaction potential in the Eulerian frame 
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Ueu (^) = ^ / V{x-y)p (f ) p [y] dy. 

Here p (x) is approximated by the ensemble average of the individual particles 
during the simulation. We also define C/La as the average collective potential 
calculated in the Lagrangian frame, L^La (^) = {U {xi))g^=g, where U (xi) is 
defined in Eq. (3) of the individual-based model. Since the rigid-body rotation 
and the single- mill state have rotational symmetry with respect to xcm, we 
evaluate Ueu and t/La after such states are reached and at position x such that 
1^— ^cm| = -R/2, where R is the swarm radius. In Fig. 12 (b), C/eu and L^La 
are shown to converge in the catastrophic regime and diverge in the H-stable 
one. In Fig. 12 (c), we investigate whether the difference between these two av- 
eraged potentials, AU = U-eu ~ Ui^a, vanishes with increasing N. Figure 12 (c) 
shows that AU indeed tends to zero for catastrophic swarms by increasing N 
but remains finite for H-stable ones. For the ensemble average to be valid in 
the H-stable regime, we may instead choose a scale that is much larger than 
the characteristic interaction lengths. Under such low resolution, the particle 
distribution can be seen as a continuous density, and Eq. (21) is then valid. 
This becomes the case of the incompressible fluids in Ref. [27], where the in- 
teraction is extremely localized, and hence, the continuum force yields a stress 
tensor as a function of the local density. However, the swarming patterns which 
we are interested in emerge on a much smaller scale. In contrast, when the 
particles are in the dispersed state, they are far away from each other; thus, 
particle-particle interaction is very weak and dominated by velocity fluctua- 
tions. The continuum force then yields a scalar pressure, which gives the gas 
dynamics equations [27]. 



5 Stability of the homogeneous solution 



That the solutions of the continuum model of Eqs. (18) and (19) relax toward 
a uniform density distribution in the H-stable regime can also be shown by 
the linear stability analysis of its homogeneous solution. Let us flrst consider 
a more general case for a 2D self-driving continuum model with a non-local 
interaction 



^-V-(p«)=0; (27) 
— {pu) + V • {puu) =f{\u\)pu-pV J V{\x-y\)p{y) dy, 

where f {\u\) is a scalar function specifying the self-driving mechanism, and 
the non-local interaction is expressed by the convolution term. For our model. 
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Fig. 12. (a) Us versus i. The upper curve represents the case of s = 10% while the 
lower curve is for s = 50%. On the upper-right corner is an illustration showing how 
the significant neighbors are defined, (b) Ueu and f/La versus i. Here a = 0.003, 
P = 0.5, Ca = 0.5, ia = 2.0, Cr = 1.0, and N = 500. (c) AU and versus N. 
i = 0.65 for the catastrophic curve while i = 0.75 for the H-stable curve. The other 
parameter values are the same as (b). 



/(|-u|) = a — P\u\ . The possible homogeneous steady state solutions can 
be written as p (x, t) = po and u (x, t) = vq v, where D is a unit vector and 
Vq can be or any of the roots of / (fo) = 0. For our Rayleigh-type dissi- 
pation, t>o = \J a/ (3. We perturb the steady state solution using p (x, t) = 
Po + 6pexp {at + iq ■ x) and u{x,t) = vqV + {Suu + Sv v) exp {at + iq ■ x), 
where 6p, 6u, Sv <^ 1 are small amplitudes. The unit vector u points to 
the direction perpendicular to v on the 2D space. The wave vector is denoted 
by q while a = a {q} represents its growth rate. By substituting this ansatz 
into Eq. (27), the dispersion relation is 



a 



^ — ipo'?sin6' —ipoqcosO 

-iqV{q)sme f {vo) 
[-iqV{q}cose f {vo) + Vof {vo) 



6u 
y6vj 



(2J 
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where a' = a + w^v ■ and V (q) is the Fourier transform of the pairwise 
interaction potential V (x) . The angle between the wave vector q and the unit 
vector u is denoted by 9. 

For the case of Vq — 0, the solution is isotropic, and we can arbitrarily choose 
the unit vector v. If the wave vector q is parallel to the arbitrarily chosen v, 
Eq. (28) reduces to 



5p 



-ipog 

/ (0) 
iqV{q) /(O) 



du 
^6vj 



and cr = / (0) or i^f (0) ± f {Oy - Apoq'^V (g)j/2. If / (0) > 0, the homoge- 
neous solution is always unstable. If / (0) < 0, the homogeneous solution is 
stable only when poq^V (q) > 0. Since po and q^ are both non-negative, the 
criterion can be reduced to 



V{q)>Q. 



(29) 



For our Raylcigh-type dissipation, / (0) = a. Since a is positive, the uniform 
density solution with zero speed is an unstable steady state solution. 

For the case oi vq satisfying / (vo) — 0, Eq. (28) becomes 



hA 



5u 
y5vj 



I 



— ipo^sin^ — ipo^cos^ 

-\qV {q) sin e 
\-iqV iq)cose Vof'{vo) J 

We thus obtain the growth rate by solving the following eigenvalue equation 



Su 



<j" - vof (vo) a'" + a'poq^V (g) - vof (vo) poq^V (g) sin^ = 0. 



./2 



(30) 



Let us consider the two cases of the wave vectors parallel and perpendic- 
ular to the v-direction. For the parallel case, i.e., ^ = 0, cr' = or u' = 

vof {vq) ± \/ {vof ("^o))^ ~ Po(fy (?) /2- On the other hand, in the perpen- 
dicular case {6 = tt/2), a' = Vof {vq) or a' = ±\J —poq'^V (q). If /' (fo) > 0, 
the homogeneous solutions are always unstable. For our Rayleigh-type dis- 
sipation, /' (vq) — —2Pvo < 0; hence, the homogeneous solution is stable 
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Fig. 13. The linear stability diagram of the swarming model of Eqs. (1)- (3). 

only when poQ^V (q) > 0, which is the same as the criterion in Eq. (29). Fur- 
ther analysis shows that for a general angle 6, Eq. (30) can be rewritten as 
- vof (vo)) {a" + poq^V (q)) + F cos^ 6 = 0, where F = v^f (vo) Poq'V (q). 
In our model, F > whenever the homogeneous solution is unstable. Thus, an 
inspection of the above equation shows that its largest root, i.e., the fastest 
growth rate, is at = 7r/2. As a result, perturbations on the direction perpen- 
dicular to the swarm velocity are the fastest growing mode, and their rate is 

-poq'^V (g) for a given q. 

Substituting Eq. (3) into Eq. (29), the linear stability criterion for our swarm- 
ing model can be exphcitly obtained as 



V (g) = 27r 



1 + g' 



./2 



3/2 



1 + Pq' 



> 0, 



(31) 



where q' = qia- Since the above criterion has to hold for all q' G M, stability 
is attained at 



ce>i if £ < 1, 
c>e if £ > 1, 
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The linear stability diagram is shown in Fig. 13. Note the close connection be- 
tween the different regimes shown here and in the H-stal:)iIity diagram of Fig. 2. 
When the homogeneous solution is linearly stable, the interaction potential is 
also H-stable. This is because the condition of Eq. (29) is also sufficient, but 
not necessary, for H-stability [29] . Further study on the dispersion relation in 
Eq. (30) reveals that cr' increases as q'^V {q) decreases, and the maximum of a' 
occurs when the minimum of q'^V (q) is reached. As a result, we are able to eval- 
uate the wavelength of the fastest growth mode and categorize the long-wave 
and the short-wave instability regions in the parameter space. Furthermore, 
we compare the fastest growth wavelength to the pattern of the fully nonlinear 
continuum model near the onset of the instability, shown in the left panel of 
Fig. 14. The simulations are initiated with a homogeneous density distribution 
and computed on a periodic domain of a 206.8 x 206.8 box. The wavenumber 
of the fastest growth mode is calculated as the minimum of Eq. (31). For the 
parameters chosen in Fig. 14, |^ = 0.121, which corresponds to a wavelength 
A = 51.87. This value matches the density aggregation patterns quite well. In 
the upper figure, a = 0; the steady state density has zero velocity, and the x-y 
directions are isotropic. In the lower figure, a ^ 0, and the velocity field of the 
swarm is initiated as m (t = 0) = ^Ja/(3y. The direction of the stripes indicates 
that the fastest growth mode is indeed perpendicular to the initial velocity, 
which is also consistent with the theoretical prediction. The results can also 
be compared to the simulation of the individual-based model by using the 
same parameter values and equivalent initial and boundary conditions. Since 
V (r) decays rapidly in r, t/ (xj) can be well approximated by including only 
the adjacent eight boxes surrounding the computational domain. The steady 
particle distributions of the individual-based simulations are shown on the 
right panel of Fig. 14. The theoretically predicted wavelength agrees with the 
patterns seen in the simulations of the continuum and the individual-based 
models. 



6 Discussion 

Soft-core interactions are widely adopted in the swarming literature [2,13,15,16,17]. 
Our investigations, with the Morse potential of Eq. (3), reveal that the com- 
monly observed core-free mill patterns only exist in the catastrophic regime 
and not in the H-stable one. In this latter case, particles arrange in rigid-body- 
like structures, similar to other H-stable interactions such as the Lennard- 
Jones potential. The morphology richness associated with soft-core catas- 
trophic potentials is one of the reasons that such interactions have been so 
broadly applied in the literature. However, a soft-core interaction, such as 
Eq. (3), does not prevent particles from occupying the same space, which is an 
unphysical situation. This can be resolved in two ways. One is that animals 
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Fig. 14. Left panel: The contours of a density distribution of the continuum model 
near the instability onset with (upper) a = and (lower) a = 1; Right panel: 
Simulations of the individual-based model using the same parameters and initial 
conditions as the left figures. The parameter values are /3 = 0.5, Ca = 0.5, Cr = 1-0, 
4 = 2.0, and 4 = 1-35. 



usually flock on a reduced dimension and thus, can use the extra dimension 
to avoid actually occupying the same space. For example, ants can crawl over 
each other; therefore, they can use z-direction to "pass through" each other 
when they flock on the x-y plane. Another way is to actually add an addi- 
tional hard-core repulsion solely to prevent overlapping. In other words, there 
is a soft-core potential that defines an equilibrium distance between particles 
and gives rise to the swarming patterns, and there is also a hard-core potential 
that specifies a forbidden distance and prevent particles from penetrating each 
other. We find that the presence of the hard-core potential affects the swarm- 
ing pattern only when N is large enough, and hence the equilibrium distance 
between nearest neighbors, determined by the collective soft-core interaction, 
collapses to the vicinity of the hard-core forbidden distance. Otherwise, flocks 
exhibit the same soft-core steady state patterns for small to moderate iV, ex- 
cept for the double-mill state, which is apparently very sensitive to hard-cores 
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and, in our simulations, are absent altogether. As increases, the equilibrium 
^NND at first decreases; the flock size increases with N only when the equilib- 
rium ^NND becomes close to the forbidden hard-core zone and cannot decrease 
further. Thus, a swarming flock at moderate N can have soft-core patterns in 
spite of the existence of a local hard-core repulsion. 

Despite the natural tendency to keep a reasonable distance between each other, 
animals may still come close and occasionally touch each other while moving 
in a biological swarm. Thus, the natural repulsive tendency can be realized as 
a soft-core repulsion while the body length of the swarming animals can be 
viewed as a hard-core forbidden zone. For biological swarms, the equilibrium 
^NND is visibly larger than the hard-core forbidden zone, which supports the 
description given in the previous paragraph. In contrast, the Lennard- Jones 
potential, used for physical systems of molecules, defines an equilibrium dis- 
tance very close to where the potential rapidly rises toward infinity. In other 
words, the equilibrium distance is nearly the same as the hard-core forbidden 
zone. Compressibility is perhaps the reason why various catastrophic patterns, 
which are not observed in the condensed phases of classical matter, can exist 
in the aggregation states of natural swarms. In artificial swarms, the hard-core 
repulsion can be understood as a collision avoidance strategy. If the distance 
to invoke the collision avoidance is much shorter than the equilibrium spacing 
between agents, various collapsing patterns shown in Ref. [28] become possible 
and might even be engineered for artificial swarming of vehicles. 



7 Summciry 

Natural swarms may switch patterns under different circumstances. While 
the changes of pattern have been understood as a result of changing individ- 
ual mobility and mutual interactions within the swarm, our individual-based 
model similarly exhibits a state transition through various swarming patterns 
by varying the parameters. The same idea can be applied to artificial swarms, 
where a group of robots can be programmed to strategically change forma- 
tions by varying the self-driving and communicating parameters of the control 
model. To analyze the stability of the emerging patterns with respect to the 
model parameters, it is advantageous to have a continuum model that can 
precisely describe the individual-based model. We illustrate a procedure to 
derive a continuum model from an individual-based model by using classical 
statistical mechanics. We show that the derived continuum model does not ap- 
proximate the individual dynamics when the interaction potential is H-stable. 
This is due to the fact that for H-stable systems, the length scale of the po- 
tential is comparable to interparticle distances, whereas in the catastrophic 
regime many particles can co-exist on a length scale comparable to the scale 
of the potential. In the catastrophic regime, the steady state solution of the 
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continuum model well matches the single-mill pattern of the individual-based 
model. The long-wave instability also shows a match to both the continuum 
and the individual-based model simulations when we theoretically analyze the 
linear stabiUty of the homogeneous solution for the continuum model. Thus, 
the continuum model may be useful for further analysis, such as the stability 
of non-trivial solutions. 
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A Derivation of the fluctuation terms 



Following Ref. [27], the momentum transport equation can be obtained by 
substituting the macroscopic momentum 

p {x, t) u (f , t) = (^^PiS {Xi - x); f^ 



into the generalized Liouville Equation, valid for non-conserved systems [42], 



dt dt 



^Y.(^- {l2Pi^ i^i - ^ij + Pk ■ Vp, (^PiS {xi - f ) j ; f ^ 
Here / is the probability density function described in Eq. (12). Since 



— ■'^xA y^Pi^ (fi - f) = — • vg^PkS (ffc - x) 

m J m 

= - Va? • ( ]6{Xk-x), 

\ m J 

/ N \ 



the transport equation can further be reduced to 
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d{pu)_j^ 



dt 



k=l 



Pkpk 

m 



6{xk- x);f) + (pkS {xk -x);f 



(A.l) 



The first term on the right hand side can be modified by noting that 



Pk 



m 



k=l 



k=l 

N N 

(Pk^ i^k - x) ; f)u + uuJ2 i^k -x);f) 

k=l k=l 
N 
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where u is the macroscopic velocity defined in Eq. (14). Eq. (A.l) then becomes 
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(A.2) 



where 



N 



We can substitute the exphcit form of ^4 from Eqs. (14) and (15) 



Pk = OiPk - P - VC/ (ffe) 



into the second term of Eq. (A.2) 
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The second term above can be further simphfied as 
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As a result, Eq. (A. 2) can be written as 



— (pit) + V • {puu) + V ■ ax = apu — 2(3Eku — 2(3qK + 2(3u ■ ax + Fy, 
which is the momentum transport equation shown in Eq. (17). 
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